Model for Nanopore Formation in Two-Dimensional Materials by Impact of Highly Charged Ions

We present a first qualitative description of the atomic dynamics in two-dimensional (2D) materials induced by the impact of slow, highly charged ions. We employ a classical molecular dynamics simulation for the motion of the target atoms combined with a Monte Carlo model for the diffusive charge transport within the layer. Depending on the velocity of charge transfer (hopping time or hole mobility) and the number of extracted electrons which, in turn, depends on the charge state of the impinging ions, we find regions of stability of the 2D structure as well as parameter combinations for which nanopore formation due to Coulomb repulsion is predicted.

I n order to exploit the special properties of two-dimensional materials for novel applications as nanoelectronic devices and sensors, it has become common to stack them into artificial layer systems of semimetals (e.g., graphene 1,2 ), semiconductors (e.g., transition-metal dichalcogenides 3−5 ), and insulators (e.g., hexagonal boron nitride 6 ) as van der Waals (vdW) heterostructures. 7, 8 Tailoring the properties of heterostructures after growth requires modification techniques with single-layer precision. Highly charged ions (HCIs) are one tool to manipulate and modify 2D heterostructures with the necessary accuracy. For example, when a freestanding vdW-MoS 2 /graphene heterostructure was recently irradiated with Xe 38+ ions, it was found that nanometer-sized pores formed only in the MoS 2 monolayer facing the ion beam while the graphene beneath remained intact. 9 The known high material selectivity in the deposition of potential energy of HCIs had been noticed earlier, 10,11 but to explore its mechanisms in more detail is essential for a targeted application in the processing of vdW heterostructures on the nanoscale: e.g., for water purification 12 or filtering CO 2 and other gases. 13,14 To test the response of 2D materials to the deposition of large amounts of potential (i.e., electronic) energy, experiments have been conducted in which HCIs were directed on such targets. 15,16 A range of materials was covered from conducting single-layer graphene (SLG), for which the ejection of numerous electrons (multiple times the initial HCI charge Q in ) and structural integrity was observed, 9,10,17 to semiconducting MoS 2 , which disintegrates around the point of ion impact in combination with the emission of only a few electrons. 18 For the latter material nanopores with a diameter of a few nanometers were observed with the mean radius depending on the initial charge state of the HCI. 11 Recently, a more systematic study on the response of freestanding fluorinated SLG (fluorographene) was performed, confirming the strong dependence of the pore radius on the electronic properties (tunable band gap, charge mobility 19−22 ) of the target. 23 It was assumed that the effect leading to pore creation could be related to bond breaking and the electrostatic repulsion of target atoms that remain charged long enough to initiate target disintegration (potential sputtering). Without reneutralization of the impact area, the repulsive Coulomb potential would quickly exceed the binding potential of the target atoms, leading to material damage that can only be avoided if charge transport in the layer proceeds fast enough. Accordingly, the high electron mobility in SLG would enable fast charge diffusion and ensure stability of the target, while MoS 2 disintegrates when excited by an HCI due to its electron mobility being about 3 orders of magnitude smaller. 24−27 To investigate this hypothesis we have set up a simulation consisting of three distinct parts which are closely intertwined. (A) The motion of atoms is modeled with a MD simulation solving Newton's equations of motion in a crystal potential and accounting for the Coulomb repulsion of positively charged ions. (B) The layer is charged due to electron transfer to the impinging projectile. (C) Positive holes left on the surface after electron transfer are distributed diffusively over the surface via stochastic hopping between neighboring lattice sites. Different charge mobilities are accounted for by the hopping time t h . While it is clear that hopping of localized charges is an inappropriate description for charge transport in conducting materials, we use the model also for SLG as a limiting case for very small hopping times: t h → 0.
Atomic units (e = m e = ℏ = 1 au) are used throughout the paper unless otherwise stated.
Molecular Dynamics Simulation. We model graphene flakes as a number of concentric rings N ring around the impact point of the HCI on the surface (middle atom, Figure 1). The distance between neighboring carbon atoms is set to the equilibrium distance at 300 K (d CC = 1.42 Å), the bond angles are 120°, i.e., a perfect honeycomb structure. Depending on the initial charge state of the projectile and the charge mobility of the target layer, the number of rings is adapted to achieve converged simulations. For example, for projectiles with a kinetic energy of 0.7 keV/amu, charge states as large as Q in = 35, and a hopping time of about 5 au we find convergence of our results starting from N ring ≥ 11 rings (∼800 atoms), which we use for most of the results presented in this letter. In our simulation atoms of the outermost ring do not move and charges hopping onto this ring are removed from the simulation and no longer contribute to Coulomb repulsion.
A variety of potentials have been used to model graphene and its thermal and mechanical properties. A popular type of potential that has been successfully used in describing 2D materials is the Stillinger−Weber (SW) potential 28 (for other possible choices see, e.g., refs 29−31). We use it in this simulation due to its small numerical cost while successfully reproducing the splitting of graphene grains into graphene nanoplatelets when strained. 32 The SW potential originally developed to model silicon 28 consists of a distance-dependent pair potential ϕ 2 and a bond potential depending on the angle between any three atoms ϕ 3 : The parameters for the potential used in this work (Table 1) are optimized to reproduce the mechanical properties of graphene. 32 Using these parameters, not only the ground-state lattice can be stabilized but also typical lattice defects such as the Stone− Wales defect (combination of rings of seven and five atoms; 33 Figure 2).
From the curvature of the potential at its minimum we calculate the oscillation frequencies in-plane and out-of-plane to determine the randomized starting positions and velocities for our MD simulations for all carbon atoms in the graphene sheet at 300 K from a microcanonical ensemble.
We solve Newton's equations of motion using a fourth-order variable time step Runge−Kutta integrator to model the motion of carbon atoms in the potential landscape defined by eq 1 and by additional unshielded positive hole charges localized on the target atoms after electron transfer from the impact area to the impinging projectile. Note that the potential (eq 1) is not altered in our simulation due to charge depletion. Furthermore, kinetic momentum transfer from the projectile onto lattice atoms is currently still neglected, as it induces mainly localized defects, e.g., a Stone−Wales defect ( Figure 2) or single vacancies, too small to be detected in typical experimental setups. 11, 34 We follow the motion of the atoms until the lattice is free of charges, i.e., charges have migrated to the "bulk" outside the outermost ring. After conclusion of the simulation, we compare the individual binding potential of each atom with its kinetic energy. Even if only a single particle's kinetic energy is larger than its binding potential, the target is counted as restructured. Again, such small changes in

Nano Letters
pubs.acs.org/NanoLett Letter the lattice structure (single vacancies) are hard to detect in experiments. Application of this criterion will therefore result in an estimate of the lower bound of the stability of SLG. Electron Transfer to Projectile. Various studies have been devoted to the simulation of electron transfer to slow HCIs in front of solids and 2D layers (e.g., refs 10 and 35−37), establishing the extremely short time scale on which the neutralization of the projectiles happens. As the positively charged ion approaches the surface, it will extract electrons from the target to excited projectile states that will eventually cascade down to less excited states in Auger or radiative decay processes. For the initial projectile charge states considered here (Q in = 30−40) a total of 90−150 electrons are estimated to be extracted from the target (sum of electrons emitted to vacuum and stabilized on the projectile 18 ).
We simulate electron transfer based on a simplified version of the classical-over-the-barrier (COB) model. 38 Accordingly, electron transfer starts at a critical distance given by with W being the work function of the material. In the current simulation we are not interested in details of the neutralization process of the HCI or the electron-emission statistics, only in the total number of electrons extracted from the impact area by the HCI. Therefore, we simplify the COB model and determine only the average distance between subsequent electron transfers Δr from the change of R c as a function of the charge state for ΔQ/Q ≪ 1. Starting with an electron transfer at R c (Q in ), we stochastically determine the distance to the next transfer as = r R ln sim c with being a uniformly distributed random number [0, 1). As the electron is captured into a highly excited state of the HCI, its screening of the nucleus is incomplete. Additionally, processes such as Auger deexcitation lead to the emission of electrons intermittently increasing the charge state of the HCI. Therefore, one transferred electron will on average lead to a reduction of the charge state ΔQ < 1. We use ΔQ as a free parameter to achieve an overall electron extraction along the trajectory of the ion as measured in experiments. 18 E.g., for Xe ions with a kinetic energy of 0.7 keV/amu and Q in = 35, an average amount of 120 electrons is extracted. 18 For this exemplary case we find the best agreement for ΔQ ≈ 0.29. This is in reasonable agreement with previous full simulations of the COB model for HCI approaching metal surfaces. 39 In our model, electrons can be transferred only from the central atom and atoms of ring 1 due to the strong distance dependence of the transfer rates 40 and each atom can supply a maximum of two electrons. Once an atom has lost both electrons, no further electron transfer to the HCI is possible unless the positive hole charges move along the layer by diffusive transport, thus reneutralizing the vicinity of the impact point and allowing for further electron transfer. With this limit to the central atoms of the flake, we are able to reproduce the small amount of extracted electrons for 2D layers with small charge mobility (MoS 2 18 ) without changing the parameter ΔQ, lending credence to the overall layout of the electron-transfer simulation.
Hopping Conduction. Transport of (localized positive) charges on the surface is accounted for by a simulation of diffusive transport by stochastic hopping of charges to neighboring lattice sites. Migration of a charge is simulated in two steps. First, a random number determines if a transfer happened during a time step Δt depending on the hopping time (average time between subsequent transfers) t h , < t t exp( / ) h , and then the new position of the charge is chosen randomly from any of the three neighboring lattice sites. If the selected site is already doubly charged, the transfer is suppressed. This part of the simulation has a single free parameter, the hopping time t h determining the diffusion constant = D d t / CC 2 h of the unbiased hopping process which, in turn, is directly proportional to the charge carrier mobility μ via the Einstein relation, eD = μk B T, with the electron (hole) charge e, the Boltzmann constant k B , and the absolute temperature T. Electron mobilities in 2D materials range from up to 200000 cm 2 /(V s) 41 for annealed and suspended single-layer graphene (typical values are 1 order of magnitude smaller) down to the mobility of monolayer MoS 2 with less than 80 cm 2 /(V s). 42 Actual velocities used in this simulation have been deduced from the band structure of the materials, accounting for the effective masses of the holes. Charge hopping continues for target atoms set in motion but still close to other target atoms, thus leading to a large fraction of neutral sputtered atoms. Currently, interactions of the hole charges with the HCI are still neglected, underestimating the diffusion speed. Furthermore, the resulting conductivity has been shown to depend not only on the charge density but also on the distribution of charges on the surface. 43 In summary, during a time step of duration Δt various processes may happen. (A) Atoms in the lattice are accelerated

Nano Letters pubs.acs.org/NanoLett
Letter due to the crystal and Coulomb potentials (B) Electron transfer from the target to the HCI occurs if at least one of the atoms in the impact region (i.e., within the first carbon ring) is less than doubly charged. (C) Localized charges on the surface migrate between atoms due to hopping conduction, affecting, in turn, A and B. We have performed a large set of calculations varying the parameters Q in and t h (hopping time or, equivalently, charge mobility of the target). For each combination of Q in and t h close to 1500 impacts were simulated.
For a fixed hopping time of t h = 7 au Figure 3 shows the final positions of carbon atoms after conclusion of the simulations for different values of Q in .
As observed in experiments, the size of the nanopores increases as a function of Q in . For this particular t h , pore formation sets in at Q in ≈ 25 where we observe restructuring of the target only in a fraction of our simulations (∼25%). For higher charge states (e.g., Q in = 40) the efficiency for pore formation approaches 100%.
The results of our simulations are summarized in a phase diagram ( Figure 4) with a relatively small transition area separating regions of stability (white area with small Q in and t h ) and pore formation (purple area with large Q in and t h ). Shaded regions indicate that both intact and altered surfaces have been found at the conclusion of the simulation. It has been shown that fluorographene becomes increasingly susceptible to pore formation with increasing fluorination. 23 The fluorination alters the band structure of SLG, introducing a band gap and reducing the charge mobility. 19−22 Increasing levels of fluorination reduce the charge mobility in fluorographene, likely aiding its susceptibilty to pore formation. At the top of Figure 4 we relate increasing hopping times to increasing degrees of fluorination (reduced mobility). While HCIinduced pores were never seen on SLG 10,17 (semimetal with very small t h ), increasing fluorination (increasing t h ) leads to pore formation with increasing probability. Irradiation of MoS 2 with its very small charge mobility (outside the range of t h covered in this work) always results in pore formation even for moderate initial charge states of the HCI. 9,11 Below a critical hopping time t c ≈ 2.5 au, surface restructuring was never observed in our simulation irrespective of the initial charge of the HCI. The transition region between stable and unstable regions loosely follows a Q t t 1/ h in c dependence. Note that the threshold found in this work provides a lower limit for the pore formation in 2D materials due to the charge-interaction effects neglected in our (unbiased) hopping simulation.
Finally, we have also analyzed the energy distribution of sputtered particles at the end of our simulation. For this analysis we concentrated on the region of the phase diagram ( Figure 4) in or close to the transition region where we expect the least influence of size effects (choice of N ring ), the chargetransport model, or the crystal potential on the final energy distribution of the particles. As can be expected due to the wide discrepancy of atom and hole velocities, most of the released carbon atoms have small kinetic energies and do not carry any charges. Once an atom is charged and starts to move due to Coulomb repulsion, its charge may still for a long time (on the time scale of the hole) be transferred to a nearby atom, leaving the starting atom without further acceleration. Furthermore, atoms are accelerated predominantly in the plane of the layer which transfers a considerable fraction of their kinetic energy on nearby atoms, thus "heating" the lattice (Figure 5a). A different heating mechanism (strong electron− phonon coupling not present in our model) has been proposed recently to explain the velocity distribution of Mo atoms released from MoS 2 supported on Au and SiO 2 . 44 In the threshold region in our study (shaded area in Figure 4) only a few atoms are sputtered (positive energies in Figure 5b) with small kinetic energies leaving behind atoms still bound in the crystal lattice (at negative energies) but with one or two neighboring atoms missing with binding energies of around −0.3 and −0.15 au, respectively. Due to the high stability of the SLG lattice, the target can be heated up locally to rather high temperatures before releasing individual atoms. Of course, this strong binding can be overcome by Coulomb forces if the hole charges are not diffusing away from the impact area quickly enough. The largest part of the energy, however, is still absorbed by the lattice and not by single sputtered particles.  To summarize, we have shown that the pore formation in clean and fluorinated graphene layers due to the impact of a highly charged ion can be qualitatively modeled with a molecular dynamics simulation for the motion of atoms in combination with Monte Carlo charge hopping for the charge conduction within the target layer. Based on the charge mobility of 2D materials, this model is able to reproduce the dependence of the pore diameter on the initial charge state of the impinging projectile or the reduction of the number of electrons extracted from materials with small charge mobility. Only for single-layer graphene, a material with high mobility, pore formation was not observed irrespective of Q in . For smaller charge mobilities disintegration of the layer (pore formation) due to repulsive forces can be expected above a threshold charge state Q in th . We find an approximate dependence on the square root of the charge mobility, Q in th . This result is stable with respect to changes in the depth of the crystal potential or the parameter governing the neutralization sequence, ΔQ (not shown).
This gives confidence that our phase diagram properly captures the qualitative dependence of pore formation in HCI−2D layer interactions on the initial charge state of the projectile and the relevant material parameter of the 2D structure, its charge mobility. Our results can therefore serve as guidelines when looking for 2D materials that are susceptible to nanopore formation by slow, highly charged ion impact.   Atoms bound in the crystal lattice have negative energies. The three distinct peaks in the distribution indicate the relative probability densities of the number of atoms that still have three, two, and one lattice neighbors. One au in time corresponds to 24.2 attoseconds, one au in energy to 27.2 eV.